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Abstract. Structural defects in ion crystals can be formed during a linear 
quench of the transverse trapping frequency across the mechanical instability from 
a linear chain to the zigzag structure. The density of defects after the sweep can 
be conveniently described by the Kibble-Zurek mechanism. In particular, the 
number of kinks in the zigzag ordering can be derived from a time-dependent 
Ginzburg-Landau equation for the order parameter, here the zigzag transverse 
size, under the assumption that the ions are continuously laser cooled. In a 
linear Paul trap the transition becomes inhomogeneous, being the charge density 
larger in the center and more rarefied at the edges. During the linear quench the 
mechanical instability is first crossed in the center of the chain, and a front, at 
which the mechanical instability is crossed during the quench, is identified which 
propagates along the chain from the center to the edges. If the velocity of this 
front is smaller than the sound velocity, the dynamics becomes adiabatic even in 
the thermodynamic limit and no defect is produced. Otherwise, the nucleation of 
kinks is reduced with respect to the case in which the charges are homogeneously 
distributed, leading to a new scaling of the density of kinks with the quenching 
rate. The analytical predictions are verified numerically by integrating the 
Langevin equations of motion of the ions, in presence of a time-dependent 
transverse confinement. We argue that the non-equilibrium dynamics of an ion 
chain in a Paul trap constitutes an ideal scenario to test the inhomogeneous 
extension of the Kibble-Zurek mechanism, which lacks experimental evidence to 
date. 
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Figure 1. Schematic representation of the "linear to zig-zag" phase transition 
in a homogeneous ion chain in two dimensions illustrating the structural phases 
involved. The high-symmetry phase corresponds to a linear chain, while 
the broken-symmetry phase is characterised by a doubly degenerate zig-zag 
configuration. 



1. Introduction 

Ion crystals in Paul or Penning traps represent a prominent example of a self organized 
system that is amenable to accurate experimental characterization and manipulation 
PQ. Ions confined by means of static and radio- frequency electromagnetic potentials 
reach crystallization when laser cooled. Crystals made up from tens to millions of ions 
have been observed both in Paul traps [2| [3j 2J [5] and in Penning traps [6 . Different 
structures can be realized by varying the particle density or the trap anisotropy. In 
Ref. [2j quasi one-dimensional structures have been experimentally characterized for 
first time. Here, the first three structures encountered when decreasing the transverse 
confinement are a linear chain, a planar zigzag structure, and a helicoidal arrangement 
of the ions about the trap axis. Phase transitions separating two structures are usually 
discontinuous. In this respect, the transition from the linear chain to the zigzag 
structure constitutes an exception. In fact, it has been demonstrated numerically [8] 
and later analytically using Landau theory [S] that the transition from a linear chain 
to the zigzag is of second order. The transition can be induced either by increasing the 
ion density or by decreasing the transverse confinement v t so that their value exceeds 
or is below, respectively, a certain critical value determining the mechanical instability. 
Recently, the transition linear-zigzag chain has been suggested as a test-bed for many- 
body quantum effects and the creation of double- well potentials jTTj , and is well suited 
for the distinction between the nucleation of defects and quasiparticles in a symmetry 
breaking scenario |12) . 

The aim of this paper is to study the out-of-equilibrium dynamics of a linear 
chain of trapped ions when v t is lowered in time from above to below the mechanical 
instability separating the linear from the two-dimensional zigzag configuration, 
illustrated in Fig. [TJ Due to the finiteness of the sound velocity, space-like separated 
regions may develop one of the two possible zigzag orderings: odd ions up and even 
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Figure 2. Ion chain in a zigzag configuration exhibiting structural defects. The 
ions are confined in a ring trap, and we report here the distribution of charges in 
two dimensions along the ring mapped to a line for clarity. The solid line joining 
the ions serves as a guide to the eye. In the example reported there are 4 defects. 
According to the Kibble-Zurek mechanism, the average size of the domains is 
given by the correlation length at the freeze-out time, £ (see text for details). 



ions down or odd ions down and even ions up as depicted in Fig. [5] These regions are 
the analogs of magnetic domains in a ferromagnetic material and the interface between 
these domains is a structural defect. The classical and quantum properties of these 
defects and their possible use for quantum information processing were studied in [T3] . 
Defects formation can be understood from simple statistical mechanics considerations. 
If the rate of change of v t is larger than the relaxation rate of the crystal, the latter 
does not have enough time to relax to the minimum energy configuration, which is 
characterized by a perfectly ordered zigzag structure, thus exhibiting no defects and 
corresponding to the state that would have been obtained after a perfect adiabatic 
transition. In principle, for an infinite system whose phase transition is described by 
Landau theory [T3] [T3], the dynamical relaxation time scale diverges at the critical 
point and therefore no matter how slowly v t changes, there will always be proliferation 
of defects. This is the celebrated Kibble-Zurek mechanism (KZM) which can generally 
account for the nucleation of topological defects in scenarios of spontaneous symmetry 
breaking [T6l [17] . Moreover, Zurek predicted the scaling of the number of defects as a 
function of the rate of passing through the phase transition [17] . The KZM scenario 
and Zurek's prediction for the scaling of the number of defects have been verified in a 
variety of systems numerically [18] and experimentally [19] . Recently, the KZ scaling 
prediction has also been extended to quantum phase transitions [20] . 

In this work we extend the results in our proposal [23] and study the production 
of kinks in the linear to zigzag transition in ion crystals by deriving a time dependent 
Ginzburg-Landau theory for the transition. With this result, we determine the scaling 
of the number of defects with the rate of change of the transverse frequency which 
follows from the KZM and compare it with numerical simulations. In a ring trap when 
the interparticle distance is uniform, one recovers the standard KZM. By contrast, in 
a linear Paul trap the density of ions is inhomogeneous, higher at the center and more 
rarefied at the edges. When decreasing in time the transverse trap frequency, the 
value of the mechanical instability is first crossed in the center of the chain, and one 
can identify a front crossing the transition from the linear to the zigzag phase which 
moves from the center towards the edges, and whose velocity depends on the rate of 
the quench. If the velocity of this front is smaller than the sound velocity, at which 
a perturbation propagates along the chain, the dynamics becomes adiabatic. This 
result holds even in the thermodynamic limit. Otherwise, the nucleation of defects is 
partially suppressed with respect to the homogeneous case and a novel scaling of the 
density of kinks with the rate of the transition emerges in the inhomogeneous extension 
of KZM [21] [22]. The whole body of experimental results in [19] has been aimed 
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at the verification of the paradigmatic KZM in a homogeneous scenario (HKZM). 
By contrast, the inhomogeneous KZM (IKZM) [3TJ [35] lacks to date experimental 
verification, and we argue in the following that ion chains in a Paul trap constitute an 
ideal system for such goal. 

2. Ginzburg-Landau equation in presence of laser cooling for the order 
parameter 

In this section we start from the theory presented in [22] IS] > and derive a Ginzburg- 
Landau equation for the order parameter of the linear-zigzag structural transition. 
The order parameter corresponds to the position offset of the ions from the trap axis. 
For a large number of ions, using a local density approximation we can approximate the 
crystal as a continuum, so that the order parameter is a field. The theory is extended 
to the case in which the crystal motion is laser cooled, and the coupling to an external 
reservoir is described using the theoretical model developed in [Jjj). Defect formation 
by quenching the control parameter, here the transverse trap frequency, across the 
structural instability is studied by using the corresponding Langevin equation for the 
order parameter. This theoretical model allows us to estimate the density of defects 
at the end of the quench. 

2.1. Ginzburg-Landau Lagrangian for the structural phase transition 

The system we consider consists of N ions of mass m, charge Q and coordinates 
r„ = (x n ,y n , z n ) which are confined in a quasi one dimensional trap with tight 
harmonic confinement at frequency v t in the plane yz. In this section we assume 
a ring trap with large radius R, so that we may impose periodic boundary conditions 
along the x-axis. The Lagrangian describing the dynamics of the ions is 

L = T — V , (1) 

where the kinetic and potential energies take the form, respectively, 

T = im£r£, (2) 

n 

V = \m^vl + 4) + ^ E F^TT, ■ (3) 

At sufficiently low temperatures and sufficiently large transverse confinement the ions 
crystallize around the stable equilibrium points of the potential V, which are aligned 
along the x-axis, ro« = (na, 0,0), with a the interparticle distance in the x-direction 
such that a = 2nR/N. 

The stability of the linear chain along the x axis requires a transverse trap 
frequency exceeding a threshold value v[ c \ which scales with the characteristic 
frequency ujq — yj Q 2 /ma 3 . At v\°' the configuration has a structural instability, such 

(c) 

that for v t < v\ the ions are organized in a planar structure with equilibrium positions 
r 0n = (na, (-l) n (6/2)cos0, (-l) n (6/2) sin0) with 6 £ [0;2tt] the angle between the 
crystal plane and the plane xy. The structure has the form of a zigzag line, which joins 
the charges along the plane, with b the transverse size of the zigzag. An appropriate 
thermodynamic limit can be defined, letting the number of ions N — > oo while 

(c) 

keeping the value v\ and the characteristic frequency ujq fixed. This corresponds 
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to keeping fixed the linear interparticle distance a or equivalently the charge density 
[2~i] . One finds that the structural instability is a second-order phase transition, 
with control field v t (or alternatively, a) and order parameter b [5]. In particular, 
the critical value of the transverse frequency in the thermodynamic limit is given 
by u\ c ' = wo\/7C(3)/2 = (2.051 .. .)ojo, C being the Riemann-zeta function. From 
the Landau theory of the structural phase transition in Ref. 9 , one can develop 
a continuum model for the field describing the order parameter applying standard 
assumptions, thereby obtaining the Ginzburg-Landau equation. We sketch below the 
assumptions made, in order to clarify the range of validity of the model which forms 
the basis of this work. 

In the linear chain, in the harmonic limit, transverse and axial modes are 
decoupled. For periodic boundary conditions the modes have a well defined quasi- 
momentum k which takes values within the Brillouin zone. For N ions, the transverse, 
normal modes at wave vector k are the real and imaginary part of the mode 



N 

IL 

where o~ n = y n , z n is the transverse displacement of the ion at the equilibrium position 
r 0n and k £ [0,7r/a]. Using decomposition as shown in Ref. [SJ, where all the 
details of the derivation are provided, the transverse potential, obtained by expanding 
Eq. ([3]) around the equilibrium positions, takes the form 

v ~ y (0) + y (2) + u (3) + y (4) (5) 

where the label indicates the order of the expansion. In particular, 

Vi2) = l E E ^(fe)(Re{^} 2 +M^} 2 ) (6) 

fe£(0,7r/a] a=y,z 



with 

/3(fc) = ^„ 4w 2£ 1 in2 (^j (7) 

j>0 3 ^ ' 

and where we omitted to write the potential term for the axial modes. The terms 
and V^ 4 ' correspond to third- and fourth-order expansion in the fluctuations around 
the equilibrium position, and contain the coupling between the axial and radial modes. 
We now make the assumption that the ions are pinned in the axial direction and can 
only oscillate in the transverse direction and discard the coupling to the axial modes. 
In this regime, it was found in Ref. [5] that — and the fourth order term reads: 

v{4) = E E ^ukzMMmMMk (8) 

fei+fe2+fe3+fc4=0 cr,T—y,z 

with 

m>0 p— 1 

In the stability regime one finds that /3(k) is minimum at fco — K / a, corresponding 

(c) 

to the so-called zigzag mode. The critical value of the transverse trap frequency v t 
is found from the relation /3(ko) = 0. When v t < the linear chain is unstable, and 
the ions order in a zigzag structure across the trap axis. 
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Hence, sufficiently close to the critical value v[ c \ the mode of the linear chain 
which first becomes unstable and determines the equilibrium positions of the zigzag 
structure (soft mode) is the transverse mode of the linear chain with the shortest 
wave length (zigzag mode), at wave vector ko = n/a. When u t is sufficiently close to 

(c) 

the critical value v\ , an effective potential can be derived for the transverse normal 
modes "J/- with wave vector k — ko — 5k, such that aSk <C 1. The effective potential is 

composed by a quadratic component, given in Eq. ((5]), with coefficient /?(fc), Eq. (JT) 
such that [TU] 

(3(~k) . . « 6 + h 2 Sk 2 , (10) 

k — ko — 5k 

where S — v\ — vf^" 1 ano - which can take negative values depending on the value of v t . 
Within the same assumptions, we can approximate: 

A(ki,k 2 ,k3,k A ) w — r %r V -, l -— = —A, (11) 

m>0 v ' 

where A = (93C(5)/32)wq/o 2 and we obtain for the fourth order term: 

y(i) = ^ A i E n,WA (12) 

fcl+fc 2 +fc 3 + fe4=0 ^ T = !/.2 

Note that the prime index in the sum over the quasimomenta refers to the fact that 
the sum is restricted to the quasi momenta kj = ko ±<5fc, with Ska <C 1. The potential 
form V in Eq. ([5]) is valid at second order in the small parameter Ska <C 1, when the 
amplitude of the transverse oscillations is much smaller than the interparticle distance 
a. 

The limit aSk <C 1 corresponds to a long-wavelength expansion made with 
respect to the zigzag mode: the phase difference of the considered transverse modes 
from the zigzag modes varies slowly from ion to ion, and a continuum description 
can be introduced. In this limit x n — > x, where x is a continuous variable, and 
(— l) n (T„ — > ^(x), where ^(x) is now the position-dependent order parameter, which 
is related to ip? through the Fourier transform 

1 f Ht 

1% = -L / — e- iSkx r(x), ° = V,z (13) 
V N J a 

while the factor 1/a in the right-hand side gives the density of states. For a large 
number of ions the discrete sum over k vectors in the potential energy becomes an 
integral according to the rule J^k ~^ I d(Sk)Na/ (2ir). The second order potential 
term becomes: 

V m = ™Y f d(Sk)d Xl dx 2 is+h 2 Sk 2 )e -Uk( Xl - X2 ) r{xi)r{x2l (14) 

2 ^-^ J zira 

and using the integral representation of the Dirac function: 

,h^k)e tSkx = 2n,S(.r) (ir.i 

we obtain: 

vi 2) = ?E || [sr{x? + h 2 (d x r(x)) 2 } ■ (i6) 
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A similar calculations leads to: 



V (* ) = ™ A j2 f ^rvfvvf- (17) 

cr,r 

Finally, we obtain the Lagrangian L — J dxC(x) where C{x) is the Lagrangian 
density and reads 

- H a {xf -A^ a (x) 2 J2^ T ( x ) 2 ]- ( 18 ) 



We now discuss individually the parameters entering Eq. (IT51) . The parameter 

5 = v 2 i-v[ c) \ (19) 

determines whether the ground state of the system is a linear chain or a zigzag, 
depending on its sign. The parameter h = woa-i/hog 2 is a velocity, and determines 
the speed with which a transverse perturbation propagates along the chain. Finally, 
the parameter A is positive and determines the value of the order parameter when 
S < as will be shown below. The Lagrangian density we derived has the form of 
a Ginzburg-Landau equation. It is valid for the modes of the linear chain close to 
instability and extends the theory presented in [S] for the specific case in which the 
coupling between axial and transverse modes can be discarded. The minimal energy 
solution of Eq. (jT5J) is constant in space and fulfills the relation: 

^[5 + 2A{^ y2 +V 2 )] =0. (20) 

Equation (l20l) always admits the solution ip" — 0, corresponding to all the ions laying 
on the x axis. This solution is clearly stable only for 5 > 0, i.e. in the linear chain 
phase. For 5 < there is an infinite family of solutions satisfying g = ±y— 6/2A, 
with g = ^ipy" 2 + ip z2 [9j. These solutions correspond to zigzag structures with the 
same amplitude but laid on different planes. 

From the Ginzburg-Landau equation one can derive the critical exponent of the 
correlation length close to the critical point and at T = 0. In particular, under the 
assumption of a small spatial and static deformation of the order parameter at a 
certain position, the field autocorrelation function at distance x decays exponentially 
as ~ exp(— |x|/£) [14], where f = h/S 1 ^ 2 is the correlation length, such that 

1/2 

(21) 



and it diverges as £ ~ (y t — i^ c ') -1 / 2 with the corresponding critical exponent is 1/2. 
The same exponent governs the behaviour of the correlation length in the zigzag phase 
close to the critical point. 



2.2. Equation of motion of the slow modes in presence of laser cooling 

We now consider the physical situation, in which the crystal motion is laser cooled. 
For the moment, we consider the discrete distribution of ion charges forming a linear 
chain (so that the effective potential for the transverse modes is essentially quadratic). 
For Doppler cooling, an effective equation for the crystal modes can be derived, which 
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has the form of a Fokker-Planck equation for the energy distribution of the normal 
modes [25) . This can be put in terms of a Langevin equation of the form 

d 2 t r k + + = e%{t), (22) 

which is here reported for the transverse modes. The damping rate rj„ tk depends 
on the corresponding mode frequency and it can take different values depending 
on the propagation direction of the cooling lasers. The scalar £^(t) represents the 
corresponding Langevin force, such that its moments fulfill the relations 

(e*(*)> =0, (23) 

{el(t)el,(t')) = 2 natk K B T6 k , k ,6(t - t')/m, (24) 

where kb is the Boltzmann constant and T is the temperature which determines the 
thermal state due to cooling. Equation (|22l) relies on the assumption that the dynamics 
of the electronic degrees of freedom is much faster than that of the motional degrees of 
freedom so they can be adiabatically eliminated from the equations for the degrees of 
freedom of the crystal normal modes, which is generally true when all ions are driven 
by the cooling laser [26] , 

In the following we focus on the situation in which the trap frequency along the 
z-axis is much larger than that along the y-axis. In this limit we drop the u-label, and 
write an equation for the transverse modes along y and close to the instability point, 
at the transition from a linear chain to a zigzag structure in the xy plane. Here, the 
parameters n k and e k are slowly varying and can be assumed to be constant. Under 
these assumptions, the Euler-Lagrange equation for if) = ip y (x) obtained from Eq. (fT8| 
reads 

d*ip - h 2 d 2 x ip + r]dti> + 6il> + 2AiJj 3 = e(t). (25) 

. 2. 

Note that h can be obtained from the group velocity s = -^===|p at the critical point 
5 = 0. Damping introduces a characteristic relaxation time scale r, with which the 
system reaches equilibrium. Sufficiently close to the transition point, rj 2 ^> S and 

r « n/S, (26) 
see |14j and the appendix. Therefore the relaxation time r diverges when the 
transverse frequency v t approaches the critical value u[ as t ~ (y t — v\ The 
same scaling is valid also in the zigzag phase 5 < 0. The quantities £ and r and their 
scaling with 6 are crucial for determining the scaling of the defect production during 
a quench of the transverse frequency as explained in the next section. 



2. 3. Quenching the trap frequency through the critical value 

Within the Ginzburg-Landau description we now assume that the transverse trap 
frequency v t undergoes a change in time in the interval [—tq,Tq], such that 

v t = ^l c)2 +S(t), (27) 

and 

S(t) = -5 — , (28) 
T Q 

with Sq > and <5o <C v[ c ^ ■ The transverse trap frequency value is swept through 
the mechanical instability of the linear-zigzag chain, such that 8{—tq) = 5q > and 
S(tq) = —Sq < 0. Correspondingly, the equation for the field now reads 

d?ip - h 2 d 2 x i> + ndti> + s(t)i> + 2A^ 3 = s(t), (29) 
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i = r(i) (30) 



where for a given time t the ground state and thermodynamic properties at the phase 
transition are well defined. The time-dependence of S(t) gives now a non-equilibrium 
problem. In principle, this should be solved considering the differential equation with 
time-dependent parameters. Nevertheless one might reach an understanding of the 
essential features by separating it into two domains. In the first domain, one can use 
the equilibrium solution provided that one may assume the adiabatic approximation, 
namely, \5(t)/6(t)\ ^> r(t), where r(t) is the relaxation time for the equilibrium 
situation at the given value 5 = 5(t). The second domain is known as the impulse 
region, where as a result of the critical slowing down the order paramer is assumed 
to be frozen, ceasing to react to the external quench [17]. This stage is separated 
from the adiabatic regime by the freeze-out time scale t, i.e., a time scale in which the 
adiabatic condition ceases to be valid. This time scale can be estimated by setting 
|<5(t)/<5(t)| = r(t), which clearly sets a lower limit to t. On this time-scale the order 
parameter is correlated over domains of characteristic length £ = £(<5(t)) given by the 
correlation length at the freeze-out time. The density of defects d (number of defects 
over the total number of ions) can then be estimated by the relation d ~ l/£. 

We now determine the density of defects for the chosen time variation of the 
parameter S, Eq. ([28]) . In this case, an estimate for the freeze-out time t is given by 
the instant of time t at which the rate of change of S equals the relaxation time: 

m 

8(i) 

Equation (I30|) allows for simple solutions in two specific limits, which will be analyzed 

in this paper. In the limit, in which \J S(i) <C i] or equivalently rf 3> Sq/tq, the 
damping overcomes the oscillations associated with the frequency at the freeze-out 

time scale, \J S(i). Note that this condition is imposed at the freeze-out-time before 
crossing the transition. It is then when the properties of the broken symmetry phase 
are determined. We denote this regime by "overdamped limit" , following the definition 
introduced in Ref. 28 . In this case, using Eqs. ([28]) and ([26]) one finds t = ^/totq 
with t = rj/8 , which sets 5(t) = ?]S / tq . The density of defects in the overdamped 
limit takes the form 

.4 = ^^)'" (3D 
to auj \tq J 

where we used Eq. (f2~Tj). assuming that the temperature is sufficiently low to neglect 
finite temperature effects. Another limit, which we will consider, is the one in which 

\J S(i) 3> ?y, namely, the Ginzburg-Landau equation is the one of an underdamped 
oscillator at the freeze-out time. We denote this regime by "underdamped" limit 
according to Ref. [28]- In this case, i = (tqTq) 1 / 3 with t = S(i) = (5 /tq) 2 / 3 , 

and the density of defects is given by 

*, > i_L(*y. (32) 

tu auj o \tqJ 

The algebraic scaling in Eq. (|52"j) coincides with the one found in Ref. [35]. The 
different power scaling of So, with respect to Ref. [35J, is due to the fact that in 
the present expressions we have written explicitly the proportionality factors (which 
in [28] are summarized in the parameter £o)- We finally note that, the fact that the 
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log[l/vx Q ] 

Figure 3. Scaling of the density of defects for in an ion chain in a ring trap as 
a function of the rate of quenching the frequency of the transverse confinement. 
The fit is logd = -0.889 + 0.239 log r (regression coefficient 0.994) for N = 50, 
and it involves average over 200 realizations. The parameters are r\ = 185a>o, 
u[ c) = 2.05cjo, S = 0.68wg, e = 2.5 X 10" 3 au^ 2 . I = Na is the length of the trap 
and e is the parameter which controls the strength of the noise which is chosen by 
adding the following term to the equation of motion eN(0, 1)V At, where At is the 

time step. This noise corresponds to a temperature of ksT = 6.28 X 10~ 6 — - — °-, 
which in our case corresponds to a temperature of the order of 20 phonons. Note 
that the ratio of the number of defects over the total number of ions N is related 
to the spatial density of defects by a factor l/N. 



density of defects exhibits a power-law behaviour as a function of tq is obviously 
due to the specific choice of the quench as a function of time and of the considered 
regime (determined by the ratio ij/8(t)). The result can be understood as the domain 
size, determining the density of defects is given by the speed at which a perturbation 
propagates along the chain, awrj, multiplied by the relaxation time scale. 

The scenario just discussed describes the dynamics of the structural phase 
transition for ions placed on a ring with equilibrium positions ro n , where the 
interparticle spacing is homogeneous. We now verify numerically the HKZM prediction 
by integrating numerically the Euler-Lagrange equations obtained by minimizing the 
Lagrangian in Eq. ([T]) for a finite number of ions, and whose motion is damped in 
presence of a Langevin force. Assuming a ring confinement and pinning of one ion, 
the relevant degrees of freedom in the critical region are the transverse coordinates 
y n . At t — the ions are in a linear chain configuration with (y n ) = 0. The system 
is then driven by a linear quench of the form given in Eq. (|28[) , and the density of 
defects d is computed at some asymptotic time, after which d remains practically 
constant. A typical evolution instance with 4 defects is shown in Fig. [2] These defects 
resemble the non-massive kinks of the Frenkel-Kontorova model with a transversal 
degree of freedom that can be described by an effective 4 theory for the translational 
displacement [29|. The classical and quantum behavior of these kinks was studied in 

In the numerics, for each set of parameters, the density of defects is averaged over 
many different realizations such as the one in Fig. [2j We then study the dependence of 
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d on tq. A least-squares fit to the list of data is used to compute the scaling exponent. 
Figure [3J shows the scaling of the density of defects d as a function of the frequency 
quench time tq in the overdamped regime, in agreement with the homogeneous KZM. 
Deviations from KZM occur for fast quenches, for which the density of defects saturates 
as a result of the interactions between different kinks. In addition, finite size effects 
limits also the validity of the KZM scaling for very slow quenches that can lead to an 
adiabatic dynamics. Nucleation of defects is expected to be suppressed whenever the 
correlation length at the freeze-out time scale exceeds the size of the system, £ > Na, 
both in the underdamped and overdamped regimes. 

3. Ion chain inside a linear Paul trap: inhomogeneous effects 

The standard scaling of topological defects for homogeneous phase transitions should 
be revised whenever the quench is local or there is a spatial dependence of the relative 
frequency S [5T] . In this section we discuss the dynamics of the structural phase 
transition for ion chains with open boundaries, for which both transverse and axial 
trapping potentials are harmonic. The ion chain in a linear Paul trap is characterized 
by a linear charge distribution which is inhomogeneous, with ion density increasing 
towards the center of the trap [50]. This leads to a stronger Coulomb repulsion for 
the ions near the center of the chain. Correspondingly, the radial short-wavelength 
modes have amplitude which is larger at the center, and vanishes at the edges p4] . 
The structural instability is hence first visible at the center of the trap, such that 
at the critical value of the transverse frequency the central ions are displaced from 
the trap axis forming a zigzag chain, as illustrated in Fig. 2J At lower values of the 
transverse frequency the number of ions which are displaced from the center increases 
towards the edges of the chain, until all the chain is in a planar configuration [7 . For 
a sufficiently long chain, one could associate with this behaviour a spatially-dependent 
critical frequency, determining the transition to the zigzag, such that it is largest at 
the center and smaller at the edges. Correspondingly, if a quench of the transverse 
frequency is applied, the transition point is crossed at different instant of times along 
the chain (from the center to the edges). The velocity, at which this front propagates, 
is the so-called front velocity 22, 23 . In this case, the ratio between the front velocity 
and the sound velocity determines nucleation of defects. In the following, we compute 
the number of kinks formed in an ion crystal in a linear Paul trap in such scenario. 
The mechanism resembles the formation of solitons in a cigar-shaped Bose-Einstein 
condensate recently discussed by Zurek [2"T] . 

More precisely, we now consider the case, in which the ions are also trapped in 
the x-direction by a harmonic potential of frequency v, such that the aspect ratio v/v t 
is sufficiently small to allow for low-dimensional crystalline structures. We will now 
consider the case in which the number N of particles is finite, but it is sufficiently 
large to allow for a local density approximation. In this limit and away from the chain 
edges the linear density n(x) is well approximated by the function |30| 



with L the half-length of the chain and x the distance from the center. Within this 
treatment, the interparticle spacing a(x) is a slowly-varying function of the position, 
such that a{x) — l/n(x). In the thermodynamic limit, in which a(0) is fixed as the 
number of particles goes to infinity, N — > oo, one recovers the statistical mechanics 




(33) 



Spontaneous nucleation of structural defects in inhomogeneous ion chains 



12 




Figure 4. Sequence of classical ground states of a harmonically trapped ID 
ion crystal for decreasing values of the transverse trap frequency vt (from top 
to bottom), across the mechanical instability from a linear to a zigzag chain. 
Due to the harmonic axial potential, the density of ions in the center is larger. 
Correspondingly, this is the first region where the zigzag structure is formed when 
decreasing vt- 



and dynamical properties of the ion chain in a infinite-radius ring trap [24[ [9] . For N 
finite, the transition from a linear to a zigzag chain can be estimated with the value 
iff fa 3Nv/ (4-^/log iV) , which was found by taking only nearest- neighbours coupling 



and where the corrections scale with powers of 1/ log N [24] . 
order unity, it corresponds to the relation 



,(c)2 



2 ma 



(0) s 



Apart from a factor of 



(34) 



with o(0) = l/n(0). 

In order to determine a Ginzburg-Landau equation for this case, we first assume 
that the spatial variation is very slow, such that sufficiently far away from the edges, 
the length scale over which the intcrparticlc distance changes is much larger than the 
coarse graining length Sx 3> a(x) with which we study the dynamics: 

'" l »a(x) (35) 



Sx = a(x)/ 



(fa- 



in this limit, we can make a slowly- varying ansatz for the short- wavelength eigenmodes 
of the linear chain, such that for a given eigenmode we extend the treatment for the 
homogeneous case starting from Eq. (JIJ and write a n = a n e na , with a n slowly- 
varying amplitude |24j . Within a local-density approximation, the ions contained in a 
region of size Sx and centered at position x, become unstable at the position-dependent 
critical transverse frequency given by 



(xf 



|C(3)— 

2 ma(x)° 



(36) 
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which has the same form as the one for the homogeneous case but with the appropriate 
interparticle distance a{x). Moreover, extending to the non- homogeneous case the 
Lagrangian in Eq. ([T8")) . we find the Lagrangian L' = J dx£'(x), with the Lagrangian 
density 

C'(x) = l -p{x) [{dMx)f ~ h{xf (d x ^{x)) 2 

- 6(x)i>(xf - A(x)ij(x) 4 ] . (37) 

In Eq. p7[) . p(x) = mn(x) is the linear mass density, and the spatial dependence of the 
coefficients 5(x), h(x), and A(x) is found by using the position-dependent interparticle 
distance a(x) in the corresponding formulas for the homogeneous case. When deriving 
the Lagrangian density of Eq. (|37|) we have neglected the coupling between axial and 
transverse degrees of freedom, assuming very small fluctuations about the critical 
point. 

In order to derive the Euler-Lagrange equation for the field, we assume that 
the linear density, and hence the interparticle distance a(x), does not depend on the 
value of the transverse trap frequency and thus remains constant when quenching 
i>t through the critical point. Terms proportional to the spatial gradient da/dx 
are also neglected, assuming that sufficiently far away from the edges the inequality 
|^"(a;)| 3> (x)a' (x) / a(x)\ holds (which is consistent with the approximation made 
for deriving the Lagrange density). Within this limit, the equation of motion for the 
field tjj = 4>( x , t) reads 

dfy - h{xfdli) + r)d t ip + 5{x, t)ip + 2A{x)^ 3 = s(t), (38) 



where now 



S(x,t) = v t {ty - v c t {x) 



= ^\0f-^\ x r-S o -, (39) 

T Q 

and we considered the time-dependence given in Eq. (|28[) . 

We now consider defect formation when the value of the transverse frequency is 
quenched through the critical point. This will happen in the central region of the chain 
first, giving rise to a propagating front along the axis, whose coordinates (xF,tp) 
satisfy S(xF,tp) = 0. The front velocity vf, at which the instability propagates, 
can be found by taking the ratio between the characteristic length of the control 
parameter, (d x S(x, t)/5(x, t))~ , over the characteristic time scale at which it changes, 
(dt5(x, t)/S(x, t))~ . It takes the form giving 

d x 5(x,t) 

An explicit dependence on the physical parameter can be found using the spatially- 
dependent critical frequency in the expression 

^^(OMl-* 2 ] 3 , (41) 

with X = x/L. For the front velocity one obtains 

dv 2 c (x) "' 



<5o 

v F ~ — 



dx 



U " 1 (1-X 2 )- 2 - (42) 



, 6^(0)V Q 1*1 

Whenever the transition is homogeneous, vp becomes infinite and the standard KZM 
applies, allowing for the nucleation of defects in the whole system. Otherwise, as 
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we show in the following, nucleation of kinks will only take place in an restricted 
fraction of the chain, with length 2X^L, where the front velocity vp is larger than the 
characteristic velocity v x with which a perturbation propagates along the chain at the 
freeze-out time. We will find that the scaling of kinks density with the quenching rate 
will be in this case different from the scaling found in the homogeneous case. 

In order to determine the density of defects, we now evaluate the velocity v x with 
which a perturbation propagates along the chain at the freeze-out time. An upper- 
bound for v x can be found by considering the ratio of the correlation length £ x and 
the relaxation time t x at t, 

v x ~k- (43) 

Tx 

where the x-subindex underlines the spatial dependence due to the inhomogeneous 
nature of the system. This also corresponds to the speed of sound at the freeze- 
out point at an energy which is one over the relaxation time. In order to compute 
it, first note that the relative frequency can be written with reference to tp as 
5(x, t) = —5o(t — tp)/rQ, where the spatial dependence is encoded in t F = t[v\ c \q) 2 — 
v\°\x) 2 ]/5q. One can find the instant t, relative to tF, at which the dynamics stops 
being adiabatic by equating the time scale 5/5 to the relaxation time t x = r//5(x,t). 
In the overdamped regime, t = (tjtq/So) 1 ^ 2 , which sets the freezed-out correlation 

length £r = auo/\J \S(x, t)| = aLOa(rj5a/TQ)^ 1,>4 . Hence, the characteristic velocity of 

a perturbation becomes v x — £, x /f x — au>$(5{)/ rfTQ) 1 / 2 . 

The key insight in the IKZM is that nucleation of defects is only expected 
whenever the front velocity is larger than v x , namely, 

v F > fix, (44) 

so that spatially separated regions of the chain are causally disconnected. The violation 
of this inequality opens the possibility of driving a truly adiabatic transition. The 
reason is that the choice of the ground state in the broken symmetry phase is not 
independent in different regions of the system whenever vp is small enough with 
respect to fi x . As a result the inhomogeneous nature of the transition allows for an 
adiabatic crossing, where a zigzag chain results after the quench with no defects. This 
holds even in thermodynamic limit and dramatically differs from the HKZM where 
defects would always be expected for infinite systems. 

In the overdamped regime, nucleation of kinks is still possible whenever 

— =Ao±{l-X 2 )- 2 >l, (45) 
v x \X\ 

with 

*, = -n± (^f- (46) 

This will generally be fulfilled in a limited region of the system, 2X*L, where the 
homogeneous KZM applies. One can estimate the effective size of this region where 
defect nucleation is possible by setting vp/v x — 1, and assuming I, « 1. Then, it 
follows that X* a ~ A which leads to the scaling law of the density of kinks 



Spon 
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log[l/vx Q ] log[l/vx Q ] 

Figure 5. Density of defects for a harmonically trapped ion chain as a function 
of the inverse of the sweeping rate (a) in the overdamped regime (r? = 100^), 
where the slope in the fit is 1.006 with regression coefficient 0.995. (b) in the 
underdamped regime (r] = 10u), where the slope in the fit is 1.384 with regression 
coefficient 0.997. The defects are only considered in the central = 30 ions, in 
order to minimize defect losses (TV = 50, 2000 realizations). The parameters are 

v[ c) {0) ~ 18l/, <5 = 36i/ 2 , e = 0.05«o^ 3/2 , with «g = Q 2 /mv 2 and v being the 
axial frequency. 



Though the absolute density of defects is reduced with respect to the homogeneous 
case (since X* < 1), the dependence on the quenching rate is enhanced. Nonetheless, 
such IKZM scaling breaks down for fast quenches due to the finite size of the chain, 
and for very slow quenches due to relevant defects losses as will be discussed below in 
more detail. 

An analogous description applies to the underdamped case, where the relaxation 
time is independent of the dissipation and diverges as t x = 1/ y/\8{x, t)|, leading to 
the freeze-out time t = (tq/Sq) 1 ^ 3 . In this time scale, the correlation length freezes, 
with respect to the quench time scale tq, at a value £ x = aw ( T Q/<5o) 1 ^ 3 leading to a 
uniform sound velocity v x = £ x /t x = aujo- Again, the transition remains adiabatic as 
long as 

f = A * (i _ X 2 )- 2 > 1, (48) 
v x \X\ 

where we introduced the parameter 

•A u 7-t . (49) 

6v ( t ) (0) 2 aw T Q 

Using the same argument employed in the overdamped regime, kinks are found in a 
region of size X* u ~ A u , so that 

— af 3 . ,50, 

To test the IKZM, we study numerically in Figure [S] the scaling of the density of 
kinks as a function of the quenching rate both in the overdamped and underdamped 
regimes. The results show a good agreement with the predictions in Eqs. (|47p an d 
(|50p . Nonetheless, there is a saturation of the density of kinks at high quenching rates 
due to the interactions between the kinks. 
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We further note that the ratio between the number of defects Af that nucleate 
according to the IKZM and HKZM obeys the relation 

^*-=2j£., (51) 
Nhkzm 

with a well define scaling with respect to the quenching rate, and which applies both 
in the underdamped and overdamped regimes. In addition, the HKZM is known to 
overestimate the number of defects by a numerical factor / of order unity, / ~ 4.5 in 
our case. 



3.1. Discussion 

In this section we discuss different mechanisms responsible for deviations from the 
IKZM scaling. Two new effects come into play with respect to the ring configuration: 
a) coupling between axial and transverse modes b) enhanced defect transport. We 
shall dwell on the implications of these two effects in KZM scaling. 

Axial-transverse mode coupling - The (inhomogeneous) trapping potential allows 
the ions to shift in the axial direction as the structural phase transition takes place. 
As a result, during the course of the transition the axial density of ions increases near 
the center of the trap, an effect which hinders the study of the IKZM, as expected 
from the derivation of the GLE in the thermodynamic limit and further suggested by 
numerical simulations. 

Dynamical losses of defects - In order to minimize the axial-transverse mode 
coupling it is desirable to drive the transition just in the center of the chain, so that 
the ions at the edges of the trap remain in the linear configuration. For the ground 
state, the amplitude of the transverse displacement of the ions increases monotonically 
from each of the edges of the chain towards the center of the trap. Hence, the 
effective Peierls-Nabarro potential [23 US] seen by a kink decays as one approaches 
the ends of the chain. As a consequence there is transport of defects, which provides 
a mechanism for their losses near the edges of the trap. Note that defect transport 
remains even if the longitudinal degrees of freedom of the ions are frozen on a lattice: 
the transverse motion suffices for its dynamics. A way of minimizing these defect 
losses is by making the inter-ion spacing homogeneous, which in turn makes the Peierls- 
Nabarro potential periodic along the chain. Different potentials have been proposed to 
achieve homogeneous inter- ion spacing |32) . For an ion chain in a linear Paul trap, the 
optimal axial trapping potential can be found by fitting the local Coulomb potential 
in an homogeneous chain, U({x n }) = £ n ^ n , \ Xn l Xn ,\ = k En^n' J^7] where a is the 
desired value of the inter- ion spacing. Under such axial confinement, the transition 
becomes then more homogeneous. Nonetheless, there is a local correction to the 

transverse critical frequency, i>r,eff( n ) = v r(t) 2 — J2 n ^n' \ x S -™ , | 3 wn i cn varies as a 
function of the position along the chain. Hence even when the ions are homogeneously 
spaced, the transition remains inhomogeneous. The underlying defect dynamics is also 
driven by the fact that pairs of defects with the same topological charge repel each 
other and attract otherwise. Scattering between kinks and anti-kinks (of positive 
and negative topological charge) can occur leading to their annihilation. Nonetheless 
defects stop seeing each other when they are separated by few ions (~ 5), which 
motivates the study of thes scaling with moderate densities. Under the conditions 
in Fig. ([5]) defects are stable on a scale ~ lCV^ 1 which suffices for its imaging in 
the laboratory. The mechanisms for defect losses are particularly relevant in the 
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underdamped regime, where higher scaling coefficients than those predicted by the 
KZM are observed and the fluctuations in the number of kinks increase. Generally, 
whenever defects losses are relevant the connection of the scaling of the density 
of defects with KZM can be questioned, since these process are non-universal and 
disregarded in the mechanism. As a result, the overdamped regime is single out as 
optimal for studies of the IKZM. 

We close by noticing that the minimal length of the chain required to verify the 
IKZM scaling is constrained by two finite-size effects: a) The existence of the adiabatic 
dynamics whenever the correlation length at the freeze out time, £ = £(t), equals the 
length of the part of the chain where defects are counted, and more fundamentally, 
where defects can nucleate Defects nucleate as long as 2X Sf L > £. b) The 

breakdown of the IKZM scaling at fast rates due to a saturation of the average density 
of defects. If one wishes to check the scaling by varying the quenching time between n 
and Tf (ideally ranging over few orders of magnitude), it should be possible to achieve 
the corresponding densities of defects (di and df) related as di — (r/ /ri) a df , where a 
is the IKZM scaling (1 in the overdamped regime, 4/3 in the underdamped regime), 
without jumping into the adiabatic dynamics. Long enough chains were created in 
various setups. In [2] a chain of up to 5 x 10 4 ions was created with an axial trapping 
frequency of 1 MHz. Ring configuration were created dynamically in |33j . In the same 
setup static structures of up to 50 ions were created at the temperature of 1 mK. 

4. Conclusions 

In this work we analyzed the formation of defects in a one dimensional Coulomb 
crystal during a frequency quench from the linear chain to the zigzag structure. We 
studied the cases of ions crystals in a linear Paul trap and in a ring-shaped trap and 
predicted the defects production rate. Our study shows that a Coulomb crystal is a 
particularly neat and controllable system where the homogeneous and inhomogeneous 
KZM can be tested. Despite the experimental work in [19] addressing the homogeneous 
KZM (HKZM), the inhomogeneous KZM (IKZM) [22l [23] lacks to date experimental 
verification, and the ion chains in a Paul trap are put forward here as an ideal system 
for such goal. 
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Appendix A. Relaxation time 

In order to compute the relaxation time, we consider the following autocorrelation 
function: 

9t(r) = (y n (t + r)y n (t)) = 1 £ e« k+k > a (* k (t + r)* k , (t)> (A.l) 

kk' 
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where for simplicity we just study the dynamics of the ions in one direction y and 
we inverted Eq. (Q} for the expansion of the particles coordinates in terms of normal 
modes. In the overdamped regime: rj ^> we can neglect the second order 

derivative in Eq. (j2"2")l whose solution becomes 

(t = **e~ + - / e" — <*->e fc (« ds (A.2) 

Substituting this solution in the definition of fft(r) Eq JA.lj) we get: 

1k b T 



k 



Zkr-1 ( Mlhlt \ 



(A.3) 



Now for fixed time t, all the components with different momentum in the 
autocorrelation function gt(r) decay exponentially with the time separation r. The 
largest time scale, corresponding to the minimum frequency S = min/ £/ 3(fc), defines 
the relaxation time: 

ro = | (A.4) 

where the subscript o stresses that this result has been obtained in the overdamped 
regime. Note also that defects are localised excitations and a correction depending on 
h can be obtained from the two-point correlation function. 
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